import numpy as np
from scipy.optimize import minimize

# R(t) parameters
a_r = 3.7741
b_r = 0.5769
c_r = -0.8062
d_r = 0.0811

# A(t) parameters
a_a = 0.6509
b_a = -0.2017
c_a = -1.1990
d_a = -0.4092

# M(t) parameters
a_m = 0.120563
b_m = 0.546877
c_m = -242.873714

# V(t) parameters
a_v = -346.3630
b_v = -0.1895
c_v = -0.9567
d_v = 0.1720

# W(N_t, C, k) parameters
alpha = 0.15
lamda = 0.1
mu = 0.08

# 安全的指数函数
def safe_exp(x):
    return np.exp(np.clip(x, -100, 100))  # 限制指数函数的输入范围

def R(N_t, C):
    return N_t * a_r * safe_exp(c_r + b_r * C) - d_r * C

def F(N_t):
    return N_t * 16 / 365    

def V(N_t, C, t):
    # 对 N_t 进行适当的缩放以避免溢出
    N_t_scaled = N_t / 1e6
    return a_v + b_v * C * safe_exp(c_v * N_t_scaled) + d_v * t    

def A(N_t, C, t):
    V_value = V(N_t, C, t)
    return a_a + b_a * safe_exp(C) + c_a * N_t + d_a * V_value

def M(N_t, t):
    return c_m + b_m * N_t + a_m * t

def W(N_t, C, k, t):
    R_value = R(N_t, C)
    F_value = F(N_t)
    A_value = A(N_t, C, t)
    M_value = M(N_t, t)
    return alpha * R_value - lamda * F_value + mu * A_value + k * alpha * M_value

# 目标函数
def objective(x, t):
    N_t, C, k = x
    return -W(N_t, C, k, t)  # 最大化 W_t

# 约束条件
def constraint1(x):
    N_t, C, k = x
    return N_t  # 确保 N_t > 0

def constraint2(x):
    N_t, C, k = x
    return C  # 确保 C > 0


visitor_mean = 991194.4444444445
visitor_std = 286175.79770371574
investment_mean = 1763810.0555555555
investment_std = 714385.1499613319

# 初始猜测值
x0 = [visitor_mean, investment_mean, 0.075]  # 选择一个在0.05和0.1之间的初始值

# 约束
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
cons = [con1, con2]

# 边界约束
bounds = [(0, None), (0, None), (0.05, 0.1)]  # N_t > 0, C > 0, 0.05 < k < 0.1

# 输入年份 t
t = 2023

# 使用拉格朗日乘子法求解
solution = minimize(objective, x0, args=(t,), method='SLSQP', constraints=cons, bounds=bounds)

# 提取结果
N_opt, C_opt, k_opt = solution.x

print(f'Optimal N(t): {N_opt * visitor_std + visitor_mean}')
print(f'Optimal C: {C_opt * investment_std + investment_mean}')
print(f'Optimal k: {k_opt}')